/*


!
!  Dalton, a molecular electronic structure program
!  Copyright (C) The Dalton Authors (see AUTHORS file for details).
!
!  This program is free software; you can redistribute it and/or
!  modify it under the terms of the GNU Lesser General Public
!  License version 2.1 as published by the Free Software Foundation.
!
!  This program is distributed in the hope that it will be useful,
!  but WITHOUT ANY WARRANTY; without even the implied warranty of
!  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
!  Lesser General Public License for more details.
!
!  If a copy of the GNU LGPL v2.1 was not distributed with this
!  code, you can obtain one at https://www.gnu.org/licenses/old-licenses/lgpl-2.1.en.html.
!

!

*/
/*-*-mode: C; c-indentation-style: "bsd"; c-basic-offset: 4; -*-*/
/* fun-ft97ax.c:

   Automatically generated code implementing FT97AX functional and
   its derivatives. It is generated by func-codegen.pl being a part of
   a "Automatic code generation framework for analytical functional
   derivative evaluation", Pawel Salek, 2005

    This functional is connected by making following changes:
    1. add "extern Functional ft97axFunctional;" to 'functionals.h'
    2. add "&ft97axFunctional," to 'functionals.c'
    3. add "fun-ft97ax.c" to 'Makefile.am', 'Makefile.in' or 'Makefile'.

    This functional has been generated from following input:
    ------ cut here -------
rho:  rhoa + rhob;
grad: sqrt(grada*grada + gradb*gradb + 2*gradab);
zeta: (rhoa-rhob)/(rhoa+rhob);

beta: 0.00293;

xa:   grada*rhoa^(-4/3);
xb:   gradb*rhob^(-4/3);

sa:   asinh(xa^2);
sb:   asinh(xb^2);

da:   1+9*beta^2*xa^2*sa^2;
db:   1+9*beta^2*xb^2*sb^2;

Exa:  -rhoa^(4/3)*beta*xa^2/(da^(1/2));
Exb:  -rhob^(4/3)*beta*xb^2/(db^(1/2));

K(rhoa,rhob,grada,gradb,gradab):=(Exa+Exb);


    ------ cut here -------
*/

 
/* strictly conform to XOPEN ANSI C standard */
#if !defined(SYS_DEC)
/* XOPEN compliance is missing on old Tru64 4.0E Alphas and pow() prototype
 * is not specified. */
#define _XOPEN_SOURCE          500
#define _XOPEN_SOURCE_EXTENDED 1
#endif
#include <math.h>
#include <stddef.h>
#include "general.h"

#define __CVERSION__
 
#include "functionals.h"
 
/* INTERFACE PART */
static integer ft97ax_isgga(void) { return 1; } /* FIXME: detect! */
static integer ft97ax_read(const char *conf_line);
static real ft97ax_energy(const FunDensProp* dp);
static void ft97ax_first(FunFirstFuncDrv *ds,   real factor,
                         const FunDensProp* dp);
static void ft97ax_second(FunSecondFuncDrv *ds, real factor,
                          const FunDensProp* dp);
static void ft97ax_third(FunThirdFuncDrv *ds,   real factor,
                         const FunDensProp* dp);
static void ft97ax_fourth(FunFourthFuncDrv *ds,   real factor,
                          const FunDensProp* dp);
 
Functional FT97axFunctional = {
  "FT97ax",       /* name */
  ft97ax_isgga,   /* gga-corrected */
   1,
  ft97ax_read,
  NULL,
  ft97ax_energy,
  ft97ax_first,
  ft97ax_second,
  ft97ax_third,
  ft97ax_fourth
};
 
/* IMPLEMENTATION PART */
static integer
ft97ax_read(const char *conf_line)
{
    fun_set_hf_weight(0);
    return 1;
}

static real
ft97ax_energy(const FunDensProp *dp)
{
    real res;
    real rhoa = dp->rhoa, rhob = dp->rhob;
    real grada = dp->grada, gradb = dp->gradb, gradab = dp->gradab;

    real t1, t2, t3, t4;

    t1 = pow(grada,2.0);
    t2 = 1/pow(rhoa,2.666666666666667);
    t3 = pow(gradb,2.0);
    t4 = 1/pow(rhob,2.666666666666667);

   /* code */
    res = -0.00293*t3/(sqrt(7.726409999999999E-5*t3*t4*pow(asinh(t3*
        t4),2.0)+1.0)*pow(rhob,1.333333333333333))-0.00293*t1/(sqrt(7.726409999999999E-5*
        t1*t2*pow(asinh(t1*t2),2.0)+1.0)*pow(rhoa,1.333333333333333));

    return res;
}

static void
ft97ax_first_helper(real rhoa, real grada, real *res)
{    real t1, t2, t3, t4, t5, t6, t7, t8, t9, t10;
    real t11;

    t1 = pow(grada,2.0);
    t2 = 1/pow(rhoa,2.666666666666667);
    t3 = asinh(t1*t2);
    t4 = pow(t3,2.0);
    t5 = sqrt(7.726409999999999E-5*t1*t4*t2+1.0);
    t6 = 1/t5;
    t7 = pow(grada,4.0);
    t8 = 1/pow(rhoa,5.333333333333333);
    t9 = 1/sqrt(t7*t8+1.0);
    t10 = 1/pow(t5,3.0);
    t11 = 1/pow(rhoa,1.333333333333333);

   /* code */
    res[0] = 0.001465*t1*t10*t11*(-4.120751999999999E-4*t3*
        t7*t9/pow(rhoa,6.333333333333333)-2.0603759999999997E-4*t1*
        t4/pow(rhoa,3.666666666666667))+0.003906666666667*t1*t6/pow(rhoa,
        2.333333333333333);
    res[1] = 0.001465*t1*t10*t11*(3.090564E-4*t3*t8*t9*pow(grada,
        3.0)+1.5452819999999998E-4*grada*t4*t2)-0.00586*grada*t6*t11;
}

static void
ft97ax_first(FunFirstFuncDrv *ds, real factor, const FunDensProp *dp)
{
    real res[2];

    ft97ax_first_helper(dp->rhoa, dp->grada, res);
   /* Final assignment */
    ds->df1000 += factor*res[0];
    ds->df0010 += factor*res[1];


    if(fabs(dp->rhoa-dp->rhob)>1e-13 ||
       fabs(dp->grada-dp->gradb)>1e-13)
        ft97ax_first_helper(dp->rhob, dp->gradb, res);
    ds->df0100 += factor*res[0];
    ds->df0001 += factor*res[1];

}

static void
ft97ax_second_helper(real rhoa, real grada, real *res)
{
    real t1, t2, t3, t4, t5, t6, t7, t8, t9, t10;
    real t11, t12, t13, t14, t15, t16, t17, t18;
    real t19, t20, t21, t22, t23;

    t1 = pow(grada,2.0);
    t2 = 1/pow(rhoa,2.666666666666667);
    t3 = asinh(t1*t2);
    t4 = pow(t3,2.0);
    t5 = sqrt(7.726409999999999E-5*t1*t4*t2+1.0);
    t6 = 1/t5;
    t7 = 1/pow(rhoa,2.333333333333333);
    t8 = pow(grada,4.0);
    t9 = 1/pow(rhoa,5.333333333333333);
    t10 = t8*t9+1.0;
    t11 = sqrt(t10);
    t12 = 1/t11;
    t13 = 1/pow(rhoa,6.333333333333333);
    t14 = 1/pow(rhoa,3.666666666666667);
    t15 = -2.0603759999999997E-4*t1*t4*t14-4.120751999999999E-4*
        t8*t3*t12*t13;
    t16 = 1/pow(t5,3.0);
    t17 = 1/pow(rhoa,1.333333333333333);
    t18 = pow(grada,3.0);
    t19 = 1.5452819999999998E-4*grada*t4*t2+3.090564E-4*t18*
        t3*t12*t9;
    t20 = 1/pow(t5,5.0);
    t21 = 1/pow(t11,3.0);
    t22 = pow(grada,6.0);
    t23 = 1/t10;

   /* code */
    res[0] = 0.001465*t1*t15*t16*t17+0.003906666666667*t1*
        t6*t7;
    res[1] = 0.001465*t1*t19*t16*t17-0.00586*grada*t6*t17;
    res[2] = 0.001465*t1*t16*t17*(-0.0010988672*t21*t3*pow(grada,
        8.0)/pow(rhoa,12.66666666666667)+0.0010988672*t22*t23/pow(rhoa,
        10.0)+0.0037086768*t12*t3*t8/pow(rhoa,7.333333333333333)+7.554711999999998E-4*
        t1*t4/pow(rhoa,4.666666666666667))-0.009115555555556*t1*t6/
        pow(rhoa,3.333333333333333)-0.003906666666667*t1*t15*t16*t7-
        0.0021975*t1*pow(t15,2.0)*t17*t20;
    res[3] = 0.001465*t1*t16*t17*(8.241503999999999E-4*t21*
        t3*pow(grada,7.0)/pow(rhoa,11.66666666666667)-8.241503999999999E-4*
        t23*pow(grada,5.0)/pow(rhoa,9.0)-4.120751999999999E-4*grada*
        t4*t14-0.0024724512*t18*t3*t12*t13)+0.007813333333333*grada*
        t6*t7-0.001953333333333*t1*t19*t16*t7-0.0021975*t1*t15*t19*
        t20*t17+0.00293*grada*t15*t16*t17;
    res[4] = 0.001465*t1*t16*t17*(-6.181127999999999E-4*t21*
        t22*t3/pow(rhoa,10.66666666666667)+6.181127999999999E-4*t23*
        t8/pow(rhoa,8.0)+0.001545282*t1*t3*t12*t9+1.5452819999999998E-4*
        t4*t2)-0.0021975*t1*t17*pow(t19,2.0)*t20-0.00586*t6*t17+0.00586*
        grada*t19*t16*t17;

}

static void
ft97ax_second(FunSecondFuncDrv *ds, real factor, const FunDensProp* dp)
{
    real res[5];
 
    ft97ax_second_helper(dp->rhoa, dp->grada, res);

    ds->df1000 += factor*res[0];
    ds->df0010 += factor*res[1];

    ds->df2000 += factor*res[2];
    ds->df1010 += factor*res[3];
    ds->df0020 += factor*res[4];


    if(fabs(dp->rhoa-dp->rhob)>1e-13 ||
       fabs(dp->grada-dp->gradb)>1e-13)
        ft97ax_second_helper(dp->rhob, dp->gradb, res);
    ds->df0100 += factor*res[0];
    ds->df0001 += factor*res[1];

    ds->df0200 += factor*res[2];
    ds->df0101 += factor*res[3];
    ds->df0002 += factor*res[4];

}

static void
ft97ax_third_helper(real rhoa, real grada, real *res)
{
    real t1, t2, t3, t4, t5, t6, t7, t8, t9, t10;
    real t11, t12, t13, t14, t15, t16, t17, t18;
    real t19, t20, t21, t22, t23, t24, t25, t26;
    real t27, t28, t29, t30, t31, t32, t33, t34;
    real t35, t36, t37, t38, t39, t40, t41, t42;
    real t43, t44, t45;

    t1 = pow(grada,2.0);
    t2 = 1/pow(rhoa,2.666666666666667);
    t3 = asinh(t1*t2);
    t4 = pow(t3,2.0);
    t5 = sqrt(7.726409999999999E-5*t1*t4*t2+1.0);
    t6 = 1/t5;
    t7 = 1/pow(rhoa,2.333333333333333);
    t8 = pow(grada,4.0);
    t9 = 1/pow(rhoa,5.333333333333333);
    t10 = t8*t9+1.0;
    t11 = sqrt(t10);
    t12 = 1/t11;
    t13 = 1/pow(rhoa,6.333333333333333);
    t14 = 1/pow(rhoa,3.666666666666667);
    t15 = -2.0603759999999997E-4*t1*t4*t14-4.120751999999999E-4*
        t8*t3*t12*t13;
    t16 = 1/pow(t5,3.0);
    t17 = 1/pow(rhoa,1.333333333333333);
    t18 = pow(grada,3.0);
    t19 = 1.5452819999999998E-4*grada*t4*t2+3.090564E-4*t18*
        t3*t12*t9;
    t20 = 1/pow(rhoa,3.333333333333333);
    t21 = pow(t15,2.0);
    t22 = 1/pow(t5,5.0);
    t23 = pow(grada,8.0);
    t24 = 1/pow(t11,3.0);
    t25 = 1/pow(rhoa,12.66666666666667);
    t26 = pow(grada,6.0);
    t27 = 1/t10;
    t28 = 1/pow(rhoa,10.0);
    t29 = 1/pow(rhoa,7.333333333333333);
    t30 = 1/pow(rhoa,4.666666666666667);
    t31 = 7.554711999999998E-4*t1*t4*t30+0.0037086768*t8*
        t3*t12*t29+0.0010988672*t26*t27*t28-0.0010988672*t23*t3*t24*
        t25;
    t32 = pow(grada,7.0);
    t33 = 1/pow(rhoa,11.66666666666667);
    t34 = pow(grada,5.0);
    t35 = 1/pow(rhoa,9.0);
    t36 = -4.120751999999999E-4*grada*t4*t14-0.0024724512*
        t18*t3*t12*t13-8.241503999999999E-4*t34*t27*t35+8.241503999999999E-4*
        t32*t3*t24*t33;
    t37 = pow(t19,2.0);
    t38 = 1/pow(rhoa,10.66666666666667);
    t39 = 1/pow(rhoa,8.0);
    t40 = 1.5452819999999998E-4*t4*t2+0.001545282*t1*t3*t12*
        t9+6.181127999999999E-4*t8*t27*t39-6.181127999999999E-4*t26*
        t3*t24*t38;
    t41 = 1/pow(t5,7.0);
    t42 = 1/pow(t11,5.0);
    t43 = pow(grada,10.0);
    t44 = 1/pow(t10,2.0);
    t45 = pow(grada,9.0);

   /* code */
    res[0] = 0.001465*t1*t15*t16*t17+0.003906666666667*t1*
        t6*t7;
    res[1] = 0.001465*t1*t19*t16*t17-0.00586*grada*t6*t17;
    res[2] = 0.001465*t1*t31*t16*t17-0.0021975*t1*t21*t22*
        t17-0.003906666666667*t1*t15*t16*t7-0.009115555555556*t1*t6*
        t20;
    res[3] = 0.00293*grada*t15*t16*t17+0.001465*t1*t36*t16*
        t17-0.0021975*t1*t15*t19*t22*t17+0.007813333333333*grada*t6*
        t7-0.001953333333333*t1*t19*t16*t7;
    res[4] = -0.00586*t6*t17+0.00586*grada*t19*t16*t17+0.001465*
        t1*t40*t16*t17-0.0021975*t1*t37*t22*t17;
    res[5] = 0.001465*t1*t16*t17*(-0.0087909376*t3*t42*pow(grada,
        12.0)/pow(rhoa,19.0)+0.0087909376*t43*t44/pow(rhoa,16.33333333333333)+
        0.023808789333333*t23*t24*t3/pow(rhoa,13.66666666666667)-0.0208784768*
        t26*t27/pow(rhoa,11.0)-0.031226142933333*t12*t3*t8/pow(rhoa,
        8.333333333333334)-0.003525532266667*t1*t4/pow(rhoa,5.666666666666667))+
        0.030385185185185*t1*t6/pow(rhoa,4.333333333333333)+0.00879*
        t1*t21*t22*t7-0.00586*t1*t31*t16*t7+0.00549375*t1*pow(t15,
        3.0)*t17*t41+0.013673333333333*t1*t15*t16*t20-0.0065925*t1*
        t31*t15*t22*t17;
    res[6] = 0.001465*t1*t16*t17*(0.0065932032*t3*t42*pow(grada,
        11.0)/pow(rhoa,18.0)-0.0065932032*t44*t45/pow(rhoa,15.33333333333333)+
        0.0015109424*grada*t4*t30+0.017856592*t18*t3*t12*t29+0.0140105568*
        t34*t27*t28-0.0162082912*t32*t3*t24*t25)+0.00586*t1*t15*t19*
        t22*t7-0.003906666666667*t1*t36*t16*t7-0.007813333333333*grada*
        t15*t16*t7-0.018231111111111*grada*t6*t20+0.004557777777778*
        t1*t19*t16*t20+0.00549375*t1*t21*t19*t41*t17-0.004395*grada*
        t21*t22*t17-0.0021975*t1*t31*t19*t22*t17-0.004395*t1*t36*t15*
        t22*t17+0.00293*grada*t31*t16*t17;
    res[7] = 0.001465*t1*t16*t17*(-0.0049449024*t3*t42*t43/
        pow(rhoa,17.0)+0.0049449024*t23*t44/pow(rhoa,14.33333333333333)-
        0.0090656544*t8*t27*t35+0.0107139552*t26*t3*t24*t33-4.120751999999999E-4*
        t4*t14-0.0090656544*t1*t3*t12*t13)+0.007813333333333*t6*t7+
        0.00293*t1*t37*t22*t7-0.001953333333333*t1*t40*t16*t7-0.007813333333333*
        grada*t19*t16*t7+0.00549375*t1*t15*t37*t41*t17-0.0021975*t1*
        t15*t40*t22*t17-0.004395*t1*t36*t19*t22*t17-0.00879*grada*
        t15*t19*t22*t17+0.00586*grada*t36*t16*t17+0.00293*t15*t16*
        t17;
    res[8] = 0.001465*t1*t16*t17*(0.0037086768*t3*t42*t45/
        pow(rhoa,16.0)-0.0037086768*t32*t44/pow(rhoa,13.33333333333333)+
        0.0037086768*grada*t3*t12*t9+0.0055630152*t18*t27*t39-0.0067992408*
        t34*t3*t24*t38)+0.00549375*t1*t17*pow(t19,3.0)*t41-0.013185*
        grada*t37*t22*t17-0.0065925*t1*t40*t19*t22*t17+0.00879*grada*
        t40*t16*t17+0.00879*t19*t16*t17;

}

static void
ft97ax_third(FunThirdFuncDrv *ds, real factor, const FunDensProp* dp)
{
    real res[9];
 
    ft97ax_third_helper(dp->rhoa, dp->grada, res);

    ds->df1000 += factor*res[0];
    ds->df0010 += factor*res[1];

    ds->df2000 += factor*res[2];
    ds->df1010 += factor*res[3];
    ds->df0020 += factor*res[4];

    ds->df3000 += factor*res[5];
    ds->df2010 += factor*res[6];
    ds->df1020 += factor*res[7];
    ds->df0030 += factor*res[8];


    if(fabs(dp->rhoa-dp->rhob)>1e-13 ||
       fabs(dp->grada-dp->gradb)>1e-13)
        ft97ax_third_helper(dp->rhob, dp->gradb, res);

    ds->df0100 += factor*res[0];
    ds->df0001 += factor*res[1];

    ds->df0200 += factor*res[2];
    ds->df0101 += factor*res[3];
    ds->df0002 += factor*res[4];

    ds->df0300 += factor*res[5];
    ds->df0201 += factor*res[6];
    ds->df0102 += factor*res[7];
    ds->df0003 += factor*res[8];

}

static void
ft97ax_fourth_helper(real rhoa, real grada, real *res)
{
    real t1, t2, t3, t4, t5, t6, t7, t8, t9, t10;
    real t11, t12, t13, t14, t15, t16, t17, t18;
    real t19, t20, t21, t22, t23, t24, t25, t26;
    real t27, t28, t29, t30, t31, t32, t33, t34;
    real t35, t36, t37, t38, t39, t40, t41, t42;
    real t43, t44, t45, t46, t47, t48, t49, t50;
    real t51, t52, t53, t54, t55, t56, t57, t58;
    real t59, t60, t61, t62, t63, t64, t65, t66;
    real t67, t68, t69, t70, t71;

    t1 = pow(grada,2.0);
    t2 = 1/pow(rhoa,2.666666666666667);
    t3 = asinh(t1*t2);
    t4 = pow(t3,2.0);
    t5 = sqrt(7.726409999999999E-5*t1*t4*t2+1.0);
    t6 = 1/t5;
    t7 = 1/pow(rhoa,2.333333333333333);
    t8 = pow(grada,4.0);
    t9 = 1/pow(rhoa,5.333333333333333);
    t10 = t8*t9+1.0;
    t11 = sqrt(t10);
    t12 = 1/t11;
    t13 = 1/pow(rhoa,6.333333333333333);
    t14 = 1/pow(rhoa,3.666666666666667);
    t15 = -2.0603759999999997E-4*t1*t4*t14-4.120751999999999E-4*
        t8*t3*t12*t13;
    t16 = 1/pow(t5,3.0);
    t17 = 1/pow(rhoa,1.333333333333333);
    t18 = pow(grada,3.0);
    t19 = 1.5452819999999998E-4*grada*t4*t2+3.090564E-4*t18*
        t3*t12*t9;
    t20 = 1/pow(rhoa,3.333333333333333);
    t21 = pow(t15,2.0);
    t22 = 1/pow(t5,5.0);
    t23 = pow(grada,8.0);
    t24 = 1/pow(t11,3.0);
    t25 = 1/pow(rhoa,12.66666666666667);
    t26 = pow(grada,6.0);
    t27 = 1/t10;
    t28 = 1/pow(rhoa,10.0);
    t29 = 1/pow(rhoa,7.333333333333333);
    t30 = 1/pow(rhoa,4.666666666666667);
    t31 = 7.554711999999998E-4*t1*t4*t30+0.0037086768*t8*
        t3*t12*t29+0.0010988672*t26*t27*t28-0.0010988672*t23*t3*t24*
        t25;
    t32 = pow(grada,7.0);
    t33 = 1/pow(rhoa,11.66666666666667);
    t34 = pow(grada,5.0);
    t35 = 1/pow(rhoa,9.0);
    t36 = -4.120751999999999E-4*grada*t4*t14-0.0024724512*
        t18*t3*t12*t13-8.241503999999999E-4*t34*t27*t35+8.241503999999999E-4*
        t32*t3*t24*t33;
    t37 = pow(t19,2.0);
    t38 = 1/pow(rhoa,10.66666666666667);
    t39 = 1/pow(rhoa,8.0);
    t40 = 1.5452819999999998E-4*t4*t2+0.001545282*t1*t3*t12*
        t9+6.181127999999999E-4*t8*t27*t39-6.181127999999999E-4*t26*
        t3*t24*t38;
    t41 = 1/pow(rhoa,4.333333333333333);
    t42 = pow(t15,3.0);
    t43 = 1/pow(t5,7.0);
    t44 = pow(grada,12.0);
    t45 = 1/pow(t11,5.0);
    t46 = 1/pow(rhoa,19.0);
    t47 = pow(grada,10.0);
    t48 = 1/pow(t10,2.0);
    t49 = 1/pow(rhoa,16.33333333333333);
    t50 = 1/pow(rhoa,13.66666666666667);
    t51 = 1/pow(rhoa,11.0);
    t52 = 1/pow(rhoa,8.333333333333334);
    t53 = 1/pow(rhoa,5.666666666666667);
    t54 = -0.003525532266667*t1*t4*t53-0.031226142933333*
        t8*t3*t12*t52-0.0208784768*t26*t27*t51+0.023808789333333*t23*
        t3*t24*t50+0.0087909376*t47*t48*t49-0.0087909376*t44*t3*t45*
        t46;
    t55 = pow(grada,11.0);
    t56 = 1/pow(rhoa,18.0);
    t57 = pow(grada,9.0);
    t58 = 1/pow(rhoa,15.33333333333333);
    t59 = 0.0015109424*grada*t4*t30+0.017856592*t18*t3*t12*
        t29+0.0140105568*t34*t27*t28-0.0162082912*t32*t3*t24*t25-0.0065932032*
        t57*t48*t58+0.0065932032*t55*t3*t45*t56;
    t60 = 1/pow(rhoa,17.0);
    t61 = 1/pow(rhoa,14.33333333333333);
    t62 = -4.120751999999999E-4*t4*t14-0.0090656544*t1*t3*
        t12*t13-0.0090656544*t8*t27*t35+0.0107139552*t26*t3*t24*t33+
        0.0049449024*t23*t48*t61-0.0049449024*t47*t3*t45*t60;
    t63 = pow(t19,3.0);
    t64 = 1/pow(rhoa,16.0);
    t65 = 1/pow(rhoa,13.33333333333333);
    t66 = 0.0037086768*grada*t3*t12*t9+0.0055630152*t18*t27*
        t39-0.0067992408*t34*t3*t24*t38-0.0037086768*t32*t48*t65+0.0037086768*
        t57*t3*t45*t64;
    t67 = 1/pow(t5,9.0);
    t68 = 1/pow(t11,7.0);
    t69 = pow(grada,14.0);
    t70 = 1/pow(t10,3.0);
    t71 = pow(grada,13.0);

   /* code */
    res[0] = 0.001465*t1*t15*t16*t17+0.003906666666667*t1*
        t6*t7;
    res[1] = 0.001465*t1*t19*t16*t17-0.00586*grada*t6*t17;
    res[2] = 0.001465*t1*t31*t16*t17-0.0021975*t1*t21*t22*
        t17-0.003906666666667*t1*t15*t16*t7-0.009115555555556*t1*t6*
        t20;
    res[3] = 0.00293*grada*t15*t16*t17+0.001465*t1*t36*t16*
        t17-0.0021975*t1*t15*t19*t22*t17+0.007813333333333*grada*t6*
        t7-0.001953333333333*t1*t19*t16*t7;
    res[4] = -0.00586*t6*t17+0.00586*grada*t19*t16*t17+0.001465*
        t1*t40*t16*t17-0.0021975*t1*t37*t22*t17;
    res[5] = 0.001465*t1*t54*t16*t17-0.0065925*t1*t31*t15*
        t22*t17+0.00549375*t1*t42*t43*t17-0.00586*t1*t31*t16*t7+0.00879*
        t1*t21*t22*t7+0.013673333333333*t1*t15*t16*t20+0.030385185185185*
        t1*t6*t41;
    res[6] = 0.00293*grada*t31*t16*t17+0.001465*t1*t59*t16*
        t17-0.0021975*t1*t31*t19*t22*t17-0.004395*grada*t21*t22*t17-
        0.004395*t1*t36*t15*t22*t17+0.00549375*t1*t21*t19*t43*t17-
        0.007813333333333*grada*t15*t16*t7-0.003906666666667*t1*t36*
        t16*t7+0.00586*t1*t15*t19*t22*t7-0.018231111111111*grada*t6*
        t20+0.004557777777778*t1*t19*t16*t20;
    res[7] = 0.00293*t15*t16*t17+0.00586*grada*t36*t16*t17+
        0.001465*t1*t62*t16*t17-0.00879*grada*t15*t19*t22*t17-0.004395*
        t1*t36*t19*t22*t17-0.0021975*t1*t15*t40*t22*t17+0.00549375*
        t1*t15*t37*t43*t17+0.007813333333333*t6*t7-0.007813333333333*
        grada*t19*t16*t7-0.001953333333333*t1*t40*t16*t7+0.00293*t1*
        t37*t22*t7;
    res[8] = 0.00879*t19*t16*t17+0.00879*grada*t40*t16*t17+
        0.001465*t1*t66*t16*t17-0.013185*grada*t37*t22*t17-0.0065925*
        t1*t40*t19*t22*t17+0.00549375*t1*t63*t43*t17;
    res[9] = 0.001465*t1*t16*t17*(-0.117212501333333*t3*t68*
        pow(grada,16.0)/pow(rhoa,25.33333333333333)+0.117212501333333*
        t69*t70/pow(rhoa,22.66666666666667)+0.357498129066667*t3*t44*
        t45/pow(rhoa,20.0)-0.318427295288889*t47*t48/pow(rhoa,17.33333333333333)-
        0.408656502044444*t23*t24*t3/pow(rhoa,14.66666666666667)+0.312932959288889*
        t26*t27/pow(rhoa,12.0)+0.279020696533333*t12*t3*t8/pow(rhoa,
        9.333333333333334)+0.019978016177778*t1*t4/pow(rhoa,6.666666666666667))-
        0.131669135802469*t1*t6*t9-0.0293*t1*t42*t43*t7+0.03516*t1*
        t31*t15*t22*t7-0.007813333333333*t1*t54*t16*t7-0.019228125*
        t1*pow(t15,4.0)*t17*t67-0.06077037037037*t1*t15*t16*t41-0.0065925*
        t1*t17*t22*pow(t31,2.0)-0.04102*t1*t21*t22*t20+0.027346666666667*
        t1*t31*t16*t20+0.0329625*t1*t31*t21*t43*t17-0.00879*t1*t54*
        t15*t22*t17;
    res[10] = 0.001465*t1*t16*t17*(0.087909376*t3*t68*pow(grada,
        15.0)/pow(rhoa,24.33333333333333)-0.087909376*t70*t71/pow(rhoa,
        21.66666666666667)-0.007051064533333*grada*t4*t53-0.1390067008*
        t18*t3*t12*t52-0.187723146666667*t34*t27*t51+0.252922600533333*
        t32*t3*t24*t50+0.219040861866667*t57*t48*t49-0.2483439872*
        t55*t3*t45*t46)-0.021975*t1*t21*t19*t43*t7+0.01758*grada*t21*
        t22*t7+0.00879*t1*t31*t19*t22*t7+0.01758*t1*t36*t15*t22*t7-
        0.00586*t1*t59*t16*t7-0.01172*grada*t31*t16*t7+0.06077037037037*
        grada*t6*t41-0.015192592592593*t1*t19*t16*t41-0.02051*t1*t15*
        t19*t22*t20+0.013673333333333*t1*t36*t16*t20+0.027346666666667*
        grada*t15*t16*t20-0.019228125*t1*t42*t19*t67*t17+0.0109875*
        grada*t42*t43*t17+0.01648125*t1*t36*t21*t43*t17+0.01648125*
        t1*t31*t15*t19*t43*t17-0.0065925*t1*t31*t36*t22*t17-0.0021975*
        t1*t54*t19*t22*t17-0.0065925*t1*t59*t15*t22*t17-0.013185*grada*
        t31*t15*t22*t17+0.00293*grada*t54*t16*t17;
    res[11] = 0.001465*t1*t16*t17*(-0.065932032*t3*t68*t69/
        pow(rhoa,23.33333333333333)+0.065932032*t44*t70/pow(rhoa,20.66666666666667)-
        0.1477976384*t23*t48*t58+0.1697749824*t47*t3*t45*t56+0.0015109424*
        t4*t30+0.0596135456*t1*t3*t12*t29+0.105765968*t8*t27*t28-0.1491712224*
        t26*t3*t24*t25)-0.01465*t1*t15*t37*t43*t7+0.00586*t1*t15*t40*
        t22*t7+0.01172*t1*t36*t19*t22*t7+0.02344*grada*t15*t19*t22*
        t7-0.003906666666667*t1*t62*t16*t7-0.015626666666667*grada*
        t36*t16*t7-0.007813333333333*t15*t16*t7-0.004395*t1*t17*t22*
        pow(t36,2.0)-0.018231111111111*t6*t20-0.006836666666667*t1*
        t37*t22*t20+0.004557777777778*t1*t40*t16*t20+0.018231111111111*
        grada*t19*t16*t20-0.019228125*t1*t21*t37*t67*t17+0.00549375*
        t1*t21*t40*t43*t17+0.00549375*t1*t31*t37*t43*t17+0.021975*
        grada*t21*t19*t43*t17+0.021975*t1*t36*t15*t19*t43*t17-0.0021975*
        t1*t31*t40*t22*t17-0.004395*t21*t22*t17-0.004395*t1*t59*t19*
        t22*t17-0.00879*grada*t31*t19*t22*t17-0.004395*t1*t62*t15*
        t22*t17-0.01758*grada*t36*t15*t22*t17+0.00586*grada*t59*t16*
        t17+0.00293*t31*t16*t17;
    res[12] = 0.001465*t1*t16*t17*(0.049449024*t3*t68*t71/
        pow(rhoa,22.33333333333333)-0.049449024*t55*t70/pow(rhoa,19.66666666666667)+
        0.0972497472*t32*t48*t61-0.1137327552*t57*t3*t45*t60-0.0543939264*
        t18*t27*t35+0.08241504*t34*t3*t24*t33-0.0197796096*grada*t3*
        t12*t13)-0.007325*t1*t63*t43*t7+0.01758*grada*t37*t22*t7+0.00879*
        t1*t40*t19*t22*t7-0.001953333333333*t1*t66*t16*t7-0.01172*
        grada*t40*t16*t7-0.01172*t19*t16*t7-0.019228125*t1*t15*t63*
        t67*t17+0.01648125*t1*t36*t37*t43*t17+0.0329625*grada*t15*
        t37*t43*t17+0.01648125*t1*t15*t40*t19*t43*t17-0.0065925*t1*
        t36*t40*t22*t17-0.013185*grada*t15*t40*t22*t17-0.0065925*t1*
        t62*t19*t22*t17-0.02637*grada*t36*t19*t22*t17-0.013185*t15*
        t19*t22*t17-0.0021975*t1*t66*t15*t22*t17+0.00879*grada*t62*
        t16*t17+0.00879*t36*t16*t17;
    res[13] = 0.001465*t1*t16*t17*(-0.037086768*t3*t44*t68/
        pow(rhoa,21.33333333333333)+0.037086768*t47*t70/pow(rhoa,18.66666666666667)+
        0.0037086768*t3*t12*t9-0.06181128*t26*t48*t65+0.074173536*
        t23*t3*t45*t64+0.0241063992*t1*t27*t39-0.0414135576*t8*t3*
        t24*t38)-0.019228125*t1*t17*pow(t19,4.0)*t67-0.0065925*t1*
        t17*t22*pow(t40,2.0)+0.04395*grada*t63*t43*t17+0.0329625*t1*
        t40*t37*t43*t17-0.02637*t37*t22*t17-0.00879*t1*t66*t19*t22*
        t17-0.05274*grada*t40*t19*t22*t17+0.01172*grada*t66*t16*t17+
        0.01758*t40*t16*t17;

}

static void
ft97ax_fourth(FunFourthFuncDrv *ds, real factor, const FunDensProp* dp)
{
    real res[14];
 
    ft97ax_fourth_helper(dp->rhoa, dp->grada, res);

    ds->df1000 += factor*res[0];
    ds->df0010 += factor*res[1];

    ds->df2000 += factor*res[2];
    ds->df1010 += factor*res[3];
    ds->df0020 += factor*res[4];

    ds->df3000 += factor*res[5];
    ds->df2010 += factor*res[6];
    ds->df1020 += factor*res[7];
    ds->df0030 += factor*res[8];

    ds->df4000 += factor*res[9];
    ds->df3010 += factor*res[10];
    ds->df2020 += factor*res[11];
    ds->df1030 += factor*res[12];
    ds->df0040 += factor*res[13];


    if(fabs(dp->rhoa-dp->rhob)>1e-13 ||
       fabs(dp->grada-dp->gradb)>1e-13)
        ft97ax_fourth_helper(dp->rhob, dp->gradb, res);

    ds->df0100 += factor*res[0];
    ds->df0001 += factor*res[1];

    ds->df0200 += factor*res[2];
    ds->df0101 += factor*res[3];
    ds->df0002 += factor*res[4];

    ds->df0300 += factor*res[5];
    ds->df0201 += factor*res[6];
    ds->df0102 += factor*res[7];
    ds->df0003 += factor*res[8];

    ds->df0400 += factor*res[9];
    ds->df0301 += factor*res[10];
    ds->df0202 += factor*res[11];
    ds->df0103 += factor*res[12];
    ds->df0004 += factor*res[13];

}
